home *** CD-ROM | disk | FTP | other *** search
/ Enter 2004 January / enter-2004-01.iso / files / maxima-5.9.0.exe / {app} / share / maxima / 5.9.0 / src / numerical / slatec / dbesi.lisp < prev    next >
Encoding:
Text File  |  2003-02-09  |  15.8 KB  |  445 lines
  1. ;;; Compiled by f2cl version 2.0 beta 2002-05-06
  2. ;;; 
  3. ;;; Options: ((:prune-labels nil) (:auto-save t) (:relaxed-array-decls t)
  4. ;;;           (:coerce-assigns :as-needed) (:array-type ':array)
  5. ;;;           (:array-slicing t) (:declare-common nil)
  6. ;;;           (:float-format double-float))
  7.  
  8. (in-package "SLATEC")
  9.  
  10.  
  11. (let ((rttpi 0.39894228040143304) (inlim 80))
  12.   (declare (type f2cl-lib:integer4 inlim) (type double-float rttpi))
  13.   (defun dbesi (x alpha kode n y nz)
  14.     (declare (type (array double-float (*)) y)
  15.              (type f2cl-lib:integer4 nz n kode)
  16.              (type double-float alpha x))
  17.     (f2cl-lib:with-array-data (y-%data% y-%offset% y)
  18.       (declare (type f2cl-lib:integer4 y-%offset%)
  19.                (type (simple-array double-float (*)) y-%data%)
  20.                (ignorable y-%offset% y-%data%))
  21.       (prog ((temp (make-array 3 :element-type 'double-float)) (ain 0.0)
  22.              (ak 0.0) (akm 0.0) (ans 0.0) (ap 0.0) (arg 0.0) (atol 0.0)
  23.              (tolln 0.0) (dfn 0.0) (dtm 0.0) (dx 0.0) (earg 0.0) (elim 0.0)
  24.              (etx 0.0) (flgik 0.0) (fn 0.0) (fnf 0.0) (fni 0.0) (fnp1 0.0)
  25.              (fnu 0.0) (gln 0.0) (ra 0.0) (s 0.0) (sx 0.0) (sxo2 0.0) (s1 0.0)
  26.              (s2 0.0) (ta 0.0) (tb 0.0) (tfn 0.0) (tm 0.0) (tol 0.0) (trx 0.0)
  27.              (t2 0.0) (xo2 0.0) (xo2l 0.0) (z 0.0) (i 0) (ialp 0) (in 0) (is 0)
  28.              (i1 0) (k 0) (kk 0) (km 0) (kt 0) (nn 0) (ns 0) (t_ 0.0))
  29.         (declare (type f2cl-lib:integer4 ns nn kt km kk k i1 is in ialp i)
  30.                  (type (array double-float (3)) temp)
  31.                  (type double-float t_ z xo2l xo2 t2 trx tol tm tfn tb ta s2 s1
  32.                   sxo2 sx s ra gln fnu fnp1 fni fnf fn flgik etx elim earg dx
  33.                   dtm dfn tolln atol arg ap ans akm ak ain))
  34.         (setf nz 0)
  35.         (setf kt 1)
  36.         (setf ra (f2cl-lib:d1mach 3))
  37.         (setf tol (max ra 1.0000000000000002e-15))
  38.         (setf i1 (f2cl-lib:int-sub (f2cl-lib:i1mach 15)))
  39.         (setf gln (f2cl-lib:d1mach 5))
  40.         (setf elim (* 2.303 (- (* i1 gln) 3.0)))
  41.         (setf i1 (f2cl-lib:int-add (f2cl-lib:i1mach 14) 1))
  42.         (setf tolln (* 2.303 gln i1))
  43.         (setf tolln (min tolln 34.5388))
  44.         (f2cl-lib:arithmetic-if (f2cl-lib:int-sub n 1)
  45.                                 (go label590)
  46.                                 (go label10)
  47.                                 (go label20))
  48.        label10
  49.         (setf kt 2)
  50.        label20
  51.         (setf nn n)
  52.         (if (or (< kode 1) (> kode 2)) (go label570))
  53.         (f2cl-lib:arithmetic-if x (go label600) (go label30) (go label80))
  54.        label30
  55.         (f2cl-lib:arithmetic-if alpha (go label580) (go label40) (go label50))
  56.        label40
  57.         (f2cl-lib:fset (f2cl-lib:fref y-%data% (1) ((1 *)) y-%offset%) 1.0)
  58.         (if (= n 1) (go end_label))
  59.         (setf i1 2)
  60.         (go label60)
  61.        label50
  62.         (setf i1 1)
  63.        label60
  64.         (f2cl-lib:fdo (i i1 (f2cl-lib:int-add i 1))
  65.                       ((> i n) nil)
  66.           (tagbody
  67.             (f2cl-lib:fset (f2cl-lib:fref y-%data% (i) ((1 *)) y-%offset%) 0.0)
  68.            label70))
  69.         (go end_label)
  70.        label80
  71.         (if (< alpha 0.0) (go label580))
  72.         (setf ialp (f2cl-lib:int alpha))
  73.         (setf fni
  74.                 (coerce
  75.                  (the f2cl-lib:integer4
  76.                       (f2cl-lib:int-sub (f2cl-lib:int-add ialp n) 1))
  77.                  'double-float))
  78.         (setf fnf (- alpha ialp))
  79.         (setf dfn (+ fni fnf))
  80.         (setf fnu dfn)
  81.         (setf in 0)
  82.         (setf xo2 (* x 0.5))
  83.         (setf sxo2 (* xo2 xo2))
  84.         (setf etx
  85.                 (coerce (the f2cl-lib:integer4 (f2cl-lib:int-sub kode 1))
  86.                         'double-float))
  87.         (setf sx (* etx x))
  88.         (if (<= sxo2 (+ fnu 1.0)) (go label90))
  89.         (if (<= x 12.0) (go label110))
  90.         (setf fn (* 0.55 fnu fnu))
  91.         (setf fn (max 17.0 fn))
  92.         (if (>= x fn) (go label430))
  93.         (setf ans (max (- 36.0 fnu) 0.0))
  94.         (setf ns (f2cl-lib:int ans))
  95.         (setf fni (+ fni ns))
  96.         (setf dfn (+ fni fnf))
  97.         (setf fn dfn)
  98.         (setf is kt)
  99.         (setf km (f2cl-lib:int-add (f2cl-lib:int-sub n 1) ns))
  100.         (if (> km 0) (setf is 3))
  101.         (go label120)
  102.        label90
  103.         (setf fn fnu)
  104.         (setf fnp1 (+ fn 1.0))
  105.         (setf xo2l (f2cl-lib:flog xo2))
  106.         (setf is kt)
  107.         (if (<= x 0.5) (go label230))
  108.         (setf ns 0)
  109.        label100
  110.         (setf fni (+ fni ns))
  111.         (setf dfn (+ fni fnf))
  112.         (setf fn dfn)
  113.         (setf fnp1 (+ fn 1.0))
  114.         (setf is kt)
  115.         (if (> (f2cl-lib:int-add (f2cl-lib:int-sub n 1) ns) 0) (setf is 3))
  116.         (go label230)
  117.        label110
  118.         (setf xo2l (f2cl-lib:flog xo2))
  119.         (setf ns (f2cl-lib:int (- sxo2 fnu)))
  120.         (go label100)
  121.        label120
  122.         (if (= kode 2) (go label130))
  123.         (if (< alpha 1.0) (go label150))
  124.         (setf z (/ x alpha))
  125.         (setf ra (f2cl-lib:fsqrt (+ 1.0 (* z z))))
  126.         (setf gln (f2cl-lib:flog (/ (+ 1.0 ra) z)))
  127.         (setf t_ (+ (* ra (- 1.0 etx)) (/ etx (+ z ra))))
  128.         (setf arg (* alpha (- t_ gln)))
  129.         (if (> arg elim) (go label610))
  130.         (if (= km 0) (go label140))
  131.        label130
  132.         (setf z (/ x fn))
  133.         (setf ra (f2cl-lib:fsqrt (+ 1.0 (* z z))))
  134.         (setf gln (f2cl-lib:flog (/ (+ 1.0 ra) z)))
  135.         (setf t_ (+ (* ra (- 1.0 etx)) (/ etx (+ z ra))))
  136.         (setf arg (* fn (- t_ gln)))
  137.        label140
  138.         (if (< arg (- elim)) (go label280))
  139.         (go label190)
  140.        label150
  141.         (if (> x elim) (go label610))
  142.         (go label130)
  143.        label160
  144.         (if (/= km 0) (go label170))
  145.         (f2cl-lib:fset (f2cl-lib:fref y-%data% (1) ((1 *)) y-%offset%)
  146.                        (f2cl-lib:fref temp (3) ((1 3))))
  147.         (go end_label)
  148.        label170
  149.         (f2cl-lib:fset (f2cl-lib:fref temp (1) ((1 3)))
  150.                        (f2cl-lib:fref temp (3) ((1 3))))
  151.         (setf in ns)
  152.         (setf kt 1)
  153.         (setf i1 0)
  154.        label180
  155.         (setf is 2)
  156.         (setf fni (- fni 1.0))
  157.         (setf dfn (+ fni fnf))
  158.         (setf fn dfn)
  159.         (if (= i1 2) (go label350))
  160.         (setf z (/ x fn))
  161.         (setf ra (f2cl-lib:fsqrt (+ 1.0 (* z z))))
  162.         (setf gln (f2cl-lib:flog (/ (+ 1.0 ra) z)))
  163.         (setf t_ (+ (* ra (- 1.0 etx)) (/ etx (+ z ra))))
  164.         (setf arg (* fn (- t_ gln)))
  165.        label190
  166.         (setf i1 (f2cl-lib:int (abs (f2cl-lib:int-sub 3 is))))
  167.         (setf i1 (max (the f2cl-lib:integer4 i1) (the f2cl-lib:integer4 1)))
  168.         (setf flgik 1.0)
  169.         (multiple-value-bind
  170.             (var-0 var-1 var-2 var-3 var-4 var-5 var-6 var-7)
  171.             (dasyik x fn kode flgik ra arg i1
  172.              (f2cl-lib:array-slice temp double-float (is) ((1 3))))
  173.           (declare (ignore var-0 var-1 var-2 var-3 var-6 var-7))
  174.           (setf ra var-4)
  175.           (setf arg var-5))
  176.         (f2cl-lib:computed-goto (label180 label350 label510) is)
  177.        label230
  178.         (setf gln (dlngam fnp1))
  179.         (setf arg (- (* fn xo2l) gln sx))
  180.         (if (< arg (- elim)) (go label300))
  181.         (setf earg (exp arg))
  182.        label240
  183.         (setf s 1.0)
  184.         (if (< x tol) (go label260))
  185.         (setf ak 3.0)
  186.         (setf t2 1.0)
  187.         (setf t_ 1.0)
  188.         (setf s1 fn)
  189.         (f2cl-lib:fdo (k 1 (f2cl-lib:int-add k 1))
  190.                       ((> k 17) nil)
  191.           (tagbody
  192.             (setf s2 (+ t2 s1))
  193.             (setf t_ (/ (* t_ sxo2) s2))
  194.             (setf s (+ s t_))
  195.             (if (< (abs t_) tol) (go label260))
  196.             (setf t2 (+ t2 ak))
  197.             (setf ak (+ ak 2.0))
  198.             (setf s1 (+ s1 fn))
  199.            label250))
  200.        label260
  201.         (f2cl-lib:fset (f2cl-lib:fref temp (is) ((1 3))) (* s earg))
  202.         (f2cl-lib:computed-goto (label270 label350 label500) is)
  203.        label270
  204.         (setf earg (/ (* earg fn) xo2))
  205.         (setf fni (- fni 1.0))
  206.         (setf dfn (+ fni fnf))
  207.         (setf fn dfn)
  208.         (setf is 2)
  209.         (go label240)
  210.        label280
  211.         (f2cl-lib:fset (f2cl-lib:fref y-%data% (nn) ((1 *)) y-%offset%) 0.0)
  212.         (setf nn (f2cl-lib:int-sub nn 1))
  213.         (setf fni (- fni 1.0))
  214.         (setf dfn (+ fni fnf))
  215.         (setf fn dfn)
  216.         (f2cl-lib:arithmetic-if (f2cl-lib:int-sub nn 1)
  217.                                 (go label340)
  218.                                 (go label290)
  219.                                 (go label130))
  220.        label290
  221.         (setf kt 2)
  222.         (setf is 2)
  223.         (go label130)
  224.        label300
  225.         (f2cl-lib:fset (f2cl-lib:fref y-%data% (nn) ((1 *)) y-%offset%) 0.0)
  226.         (setf nn (f2cl-lib:int-sub nn 1))
  227.         (setf fnp1 fn)
  228.         (setf fni (- fni 1.0))
  229.         (setf dfn (+ fni fnf))
  230.         (setf fn dfn)
  231.         (f2cl-lib:arithmetic-if (f2cl-lib:int-sub nn 1)
  232.                                 (go label340)
  233.                                 (go label310)
  234.                                 (go label320))
  235.        label310
  236.         (setf kt 2)
  237.         (setf is 2)
  238.        label320
  239.         (if (<= sxo2 fnp1) (go label330))
  240.         (go label130)
  241.        label330
  242.         (setf arg (+ (- arg xo2l) (f2cl-lib:flog fnp1)))
  243.         (if (< arg (- elim)) (go label300))
  244.         (go label230)
  245.        label340
  246.         (setf nz (f2cl-lib:int-sub n nn))
  247.         (go end_label)
  248.        label350
  249.         (setf nz (f2cl-lib:int-sub n nn))
  250.        label360
  251.         (if (= kt 2) (go label420))
  252.         (setf s1 (f2cl-lib:fref temp (1) ((1 3))))
  253.         (setf s2 (f2cl-lib:fref temp (2) ((1 3))))
  254.         (setf trx (/ 2.0 x))
  255.         (setf dtm fni)
  256.         (setf tm (* (+ dtm fnf) trx))
  257.         (if (= in 0) (go label390))
  258.         (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
  259.                       ((> i in) nil)
  260.           (tagbody
  261.             (setf s s2)
  262.             (setf s2 (+ (* tm s2) s1))
  263.             (setf s1 s)
  264.             (setf dtm (- dtm 1.0))
  265.             (setf tm (* (+ dtm fnf) trx))
  266.            label380))
  267.         (f2cl-lib:fset (f2cl-lib:fref y-%data% (nn) ((1 *)) y-%offset%) s1)
  268.         (if (= nn 1) (go end_label))
  269.         (f2cl-lib:fset
  270.          (f2cl-lib:fref y-%data% ((f2cl-lib:int-sub nn 1)) ((1 *)) y-%offset%)
  271.          s2)
  272.         (if (= nn 2) (go end_label))
  273.         (go label400)
  274.        label390
  275.         (f2cl-lib:fset (f2cl-lib:fref y-%data% (nn) ((1 *)) y-%offset%) s1)
  276.         (f2cl-lib:fset
  277.          (f2cl-lib:fref y-%data% ((f2cl-lib:int-sub nn 1)) ((1 *)) y-%offset%)
  278.          s2)
  279.         (if (= nn 2) (go end_label))
  280.        label400
  281.         (setf k (f2cl-lib:int-add nn 1))
  282.         (f2cl-lib:fdo (i 3 (f2cl-lib:int-add i 1))
  283.                       ((> i nn) nil)
  284.           (tagbody
  285.             (setf k (f2cl-lib:int-sub k 1))
  286.             (f2cl-lib:fset
  287.              (f2cl-lib:fref y-%data%
  288.                             ((f2cl-lib:int-sub k 2))
  289.                             ((1 *))
  290.                             y-%offset%)
  291.              (+
  292.               (* tm
  293.                  (f2cl-lib:fref y-%data%
  294.                                 ((f2cl-lib:int-sub k 1))
  295.                                 ((1 *))
  296.                                 y-%offset%))
  297.               (f2cl-lib:fref y-%data% (k) ((1 *)) y-%offset%)))
  298.             (setf dtm (- dtm 1.0))
  299.             (setf tm (* (+ dtm fnf) trx))
  300.            label410))
  301.         (go end_label)
  302.        label420
  303.         (f2cl-lib:fset (f2cl-lib:fref y-%data% (1) ((1 *)) y-%offset%)
  304.                        (f2cl-lib:fref temp (2) ((1 3))))
  305.         (go end_label)
  306.        label430
  307.         (setf earg (/ rttpi (f2cl-lib:fsqrt x)))
  308.         (if (= kode 2) (go label440))
  309.         (if (> x elim) (go label610))
  310.         (setf earg (* earg (exp x)))
  311.        label440
  312.         (setf etx (* 8.0 x))
  313.         (setf is kt)
  314.         (setf in 0)
  315.         (setf fn fnu)
  316.        label450
  317.         (setf dx (+ fni fni))
  318.         (setf tm 0.0)
  319.         (if (and (= fni 0.0) (< (abs fnf) tol)) (go label460))
  320.         (setf tm (* 4.0 fnf (+ fni fni fnf)))
  321.        label460
  322.         (setf dtm (* dx dx))
  323.         (setf s1 etx)
  324.         (setf trx (- dtm 1.0))
  325.         (setf dx (/ (- (+ trx tm)) etx))
  326.         (setf t_ dx)
  327.         (setf s (+ 1.0 dx))
  328.         (setf atol (* tol (abs s)))
  329.         (setf s2 1.0)
  330.         (setf ak 8.0)
  331.         (f2cl-lib:fdo (k 1 (f2cl-lib:int-add k 1))
  332.                       ((> k 25) nil)
  333.           (tagbody
  334.             (setf s1 (+ s1 etx))
  335.             (setf s2 (+ s2 ak))
  336.             (setf dx (- dtm s2))
  337.             (setf ap (+ dx tm))
  338.             (setf t_ (/ (* (- t_) ap) s1))
  339.             (setf s (+ s t_))
  340.             (if (<= (abs t_) atol) (go label480))
  341.             (setf ak (+ ak 8.0))
  342.            label470))
  343.        label480
  344.         (f2cl-lib:fset (f2cl-lib:fref temp (is) ((1 3))) (* s earg))
  345.         (if (= is 2) (go label360))
  346.         (setf is 2)
  347.         (setf fni (- fni 1.0))
  348.         (setf dfn (+ fni fnf))
  349.         (setf fn dfn)
  350.         (go label450)
  351.        label500
  352.         (setf akm (max (- 3.0 fn) 0.0))
  353.         (setf km (f2cl-lib:int akm))
  354.         (setf tfn (+ fn km))
  355.         (setf ta
  356.                 (/ (+ (- (+ gln tfn) 0.9189385332) (/ -0.0833333333 tfn))
  357.                    (+ tfn 0.5)))
  358.         (setf ta (- xo2l ta))
  359.         (setf tb (/ (- (+ 1.0 (/ (* -1 1.0) tfn))) tfn))
  360.         (setf ain
  361.                 (+ (/ tolln (- (f2cl-lib:fsqrt (- (* ta ta) (* tolln tb))) ta))
  362.                    1.5))
  363.         (setf in (f2cl-lib:int ain))
  364.         (setf in (f2cl-lib:int-add in km))
  365.         (go label520)
  366.        label510
  367.         (setf t_ (/ 1.0 (* fn ra)))
  368.         (setf ain
  369.                 (+
  370.                  (/ tolln
  371.                     (+ gln (f2cl-lib:fsqrt (+ (* gln gln) (* t_ tolln)))))
  372.                  1.5))
  373.         (setf in (f2cl-lib:int ain))
  374.         (if (> in inlim) (go label160))
  375.        label520
  376.         (setf trx (/ 2.0 x))
  377.         (setf dtm (+ fni in))
  378.         (setf tm (* (+ dtm fnf) trx))
  379.         (setf ta 0.0)
  380.         (setf tb tol)
  381.         (setf kk 1)
  382.        label530
  383.         (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
  384.                       ((> i in) nil)
  385.           (tagbody
  386.             (setf s tb)
  387.             (setf tb (+ (* tm tb) ta))
  388.             (setf ta s)
  389.             (setf dtm (- dtm 1.0))
  390.             (setf tm (* (+ dtm fnf) trx))
  391.            label540))
  392.         (if (/= kk 1) (go label550))
  393.         (setf ta (* (/ ta tb) (f2cl-lib:fref temp (3) ((1 3)))))
  394.         (setf tb (f2cl-lib:fref temp (3) ((1 3))))
  395.         (setf kk 2)
  396.         (setf in ns)
  397.         (if (/= ns 0) (go label530))
  398.        label550
  399.         (f2cl-lib:fset (f2cl-lib:fref y-%data% (nn) ((1 *)) y-%offset%) tb)
  400.         (setf nz (f2cl-lib:int-sub n nn))
  401.         (if (= nn 1) (go end_label))
  402.         (setf tb (+ (* tm tb) ta))
  403.         (setf k (f2cl-lib:int-sub nn 1))
  404.         (f2cl-lib:fset (f2cl-lib:fref y-%data% (k) ((1 *)) y-%offset%) tb)
  405.         (if (= nn 2) (go end_label))
  406.         (setf dtm (- dtm 1.0))
  407.         (setf tm (* (+ dtm fnf) trx))
  408.         (setf km (f2cl-lib:int-sub k 1))
  409.         (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
  410.                       ((> i km) nil)
  411.           (tagbody
  412.             (f2cl-lib:fset
  413.              (f2cl-lib:fref y-%data%
  414.                             ((f2cl-lib:int-sub k 1))
  415.                             ((1 *))
  416.                             y-%offset%)
  417.              (+ (* tm (f2cl-lib:fref y-%data% (k) ((1 *)) y-%offset%))
  418.                 (f2cl-lib:fref y-%data%
  419.                                ((f2cl-lib:int-add k 1))
  420.                                ((1 *))
  421.                                y-%offset%)))
  422.             (setf dtm (- dtm 1.0))
  423.             (setf tm (* (+ dtm fnf) trx))
  424.             (setf k (f2cl-lib:int-sub k 1))
  425.            label560))
  426.         (go end_label)
  427.        label570
  428.         (xermsg "SLATEC" "DBESI" "SCALING OPTION, KODE, NOT 1 OR 2." 2 1)
  429.         (go end_label)
  430.        label580
  431.         (xermsg "SLATEC" "DBESI" "ORDER, ALPHA, LESS THAN ZERO." 2 1)
  432.         (go end_label)
  433.        label590
  434.         (xermsg "SLATEC" "DBESI" "N LESS THAN ONE." 2 1)
  435.         (go end_label)
  436.        label600
  437.         (xermsg "SLATEC" "DBESI" "X LESS THAN ZERO." 2 1)
  438.         (go end_label)
  439.        label610
  440.         (xermsg "SLATEC" "DBESI" "OVERFLOW, X TOO LARGE FOR KODE = 1." 6 1)
  441.         (go end_label)
  442.        end_label
  443.         (return (values nil nil nil nil nil nz))))))
  444.  
  445.